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ABSTRACT 

We address the dual challenge of estimating deviations from Gaussianity arising in models 
of the Early Universe, whilst retaining information necessary to assess whether a detection 
of non-Gaussianity is primordial. We do this by constructing a new statistic, the bispectrum- 
related power spectrum, which is constructed from a map of the Cosmic Microwave Back- 
ground. The estimator is optimised for primordial non-Gaussianity detection, but can also be 
useful in distinguishing primordial non-Gaussianity from secondary non-Gaussianity, such as 
may arise from unsubtracted point sources, or residuals from component separation. Extend- 
ing earlier studies we present unbiased non-Gaussianity estimators optimised for partial sky 
coverage and inhomogeneous noise associated with realistic scan strategies, but which retain 
the ability to assess foreground contamination. 

Key words: Cosmology: theory - cosmic microwave-background - large-scale structure of 
Universe Methods: analytical - Methods: statistical -Methods: numerical 



1 INTRODUCTION 

The statistical properties of fluctuations in the cosmic microwave background (CMB) radiation can be used to probe the very earliest stages 
of the Universe's history, and provide valuable information on the mechanisms which ultimately gave rise to the existence of structure within 
the Universe. This may include evidence for inflation, the process by which the rapid expansion of the Universe is thought to have arisen. 
In standard inflationary models, the fields in the early Universe should be very close to random Gaussian fields, so a detection of large 
non-Gaussianity would be highly significant, and may indicate a different history, such as warm inflation or multiple-field inflation, or a 
completely different mechanism such as those arising from topological defects. A difficulty for methods designed to detect non-Gaussianity 
in the CMB is that other processe s can contribute, such as g r avitational lens i ng, unsubtracted point sources, and imperfect subtraction of 
alact ic foreground emission (e.g. lGoldberg & Spergell Jl999h ; ICooravl j200ti) ; IVerde,& Spergell j2002h : ICastrol ( 120041) : iBabich & Pierpaolil 
2008)). The challenge therefore is to provide evidence that any detection of non-Gaussianity is primordial in origin, and not a result of 
these other effects. The aim of this paper is to provide an optimised framework not just for detecting non-Gaussianity, but for assessing the 

contributions from various sources. 

Non-Gaussianity from simplest inflationary models based on a single slowly-rolling sca lar field is typically very small llSalopek & Bond 
Il990t Il99ll : iFalk et al.|[l99l iGangui et alj|l994l ; lAcquaviva et alj|2003l ; lMaldacendl2003l) . (see iBartolo Matarrese & Riottd d2006h an d 
references there in for more details). Variants of simple inflationary models such as multiple scalar fields (Lyth, Ungarelli & Wands 2003), 
features in the inflatio nary potential, non-adiabatic fluctuations, non-standard kinetic terms, warm inflation (Gu pta. Berera & Heavensl i 2002: 
Moss & Xiongi r2007). or deviations from Bunch-Davies v acuum can all l ead to much higher level of non-Gaussianity. Early observational 
work on the bi spectrum from COBE llKomatsu et alj2002h and MAXIMA ^Santos et al.l2003l) was followed by much more accurate analysis 
with WMAP iKomatsu et alj|2003l : ICreminelli et alj|2007l ; ISpergel et alj|2007t) . With the recent claim of a detection of non-Gaussianity 
dYadav &~W andelt 200^) in the Wilkinson Microwave Anisotropy Probe 5-year (WMAP5) sky maps, interest in non-Gaussianity has obtained 
a tremendous boost. Much of the interest in primordial non-Gaussianity has focussed on a phenomenological 'lo cal /nl' parametrisation in 
terms of the perturbative non-linear coupling in the primordial curvature perturbation (Komatsu & Spergel 2 00 ll) : 

$(x) = $l(x) + fi,L(*l(x) - (*i(a:)», (1) 

where <&l(x) denotes the linear Gaussian part of the Bardeen curvature and /ml is the non-linear coupling parameter. A number of 
models have non-Gaussianity which can be approximated by this form. The leading order non-Gaussianity therefore is at the level 
of the bispectrum, or in configuration space at the three-point level. Many studies i nvolving prim ordial non-Gaussianity have used 
the bispectrum, motivated by the fact that it conta i ns all the information about Jnl i Babich 2005). It has been extensively studied 
teomatsu. Spergel & Wandelfc005l;ICreminellill2003l ; ICreminelli etai]|2006l : iMedeiros & Contaldcj|200e[ : ICabella et al. 12004 lliguori et all 
2007; Smith, S enatore & Z aldarriaga 2009J), wit h most o f these measurements providing convolved esti mates of the bispectrum. Op- 
timised 3-point estimators were introduced by iHeavensI fl993), and have been successively deve l oped iKomatsu. Spergel & Wa ndell 
120051 : ICreminelli etafl 120061 : ICreminelli. Senat ore, & Za ldarriagal 120071 : ISmith. Zahn & Dorefl200ol : ISmith & Zaldarriagd l2006t) to the 



2 Munshi & Heavens 



point where an estimator for Jnl which saturates the Cramer-Rao bound exists for partial sky coverage and inhomogeneous noise 
( Smith, Senatore & Zaldarriag3l2 009). Approximat e forms also exist for equilateral non-Gaussianity, which m ay arise in models with non- 
minimal Lagrangian with higher-derivative terms dChen. H uang & Kachru 200fi| IChen. Easther & Limll2007h . In these models, the largest 
signal comes from spherical harmonic modes with l\ ~ i% — £3, whereas for the local model, the signal is highest when one I is much 
smaller than the other two - the so-called squeezed configuration. 

Reducing the CMB data to a single lossless estimate for }nl is extremely elegant, but it suffers from the disadvantage that a single 
number loses the ability to determine the extent to which the estimate has been contaminated by non-primordial signals. What we seek to to 
perform a less aggressive data compression of the CMB data, not to a single number, but to a function, which has known expected form for 
primordial models, and for which the contributions from other signals can be estimated. The purpose of this is to be able to demonstrate that a 
non-Gaussian signal is indeed primordial, or alternatively accounted for by non-primordial signals. We do this in a way which is still optimal 
for local or equilateral /jvx models, although the formalism is general. The function we choose is the integrated cross-power spectrum of 
pair of maps constructed from the CMB data. Mathematically, it is closely related to previous estimators, but the interpretation of the output 
is different, and offers very significant advantages. 

This paper is organised as follows: section ^provides a small review of available models of primordial non-Gaussianity. Section ^3] 
relates the projected-bispectrum to the corresponding primordial bispectrum. Section ^presents the optimised bispectrum-related power 
spectrum estimator for the idealised case of an all-sky survey with homogeneous noise. This is not optimal for partial sky cover age or 
inhomogeneous noise, but is straightforward and shows the connection with the /jv l estimator of Ko matsu. Spergel & Wand elt ( 2005). Here 
we present the theoretical expectation for the local }nl model, and show the link between the Jnl estimator based on one-point statistics 
and its two-point counterpart. Section ^provides a method which can handle partial sky coverage and non-uniform noise in an approximate 
way using Monte-Carlo simulations. In section ^6]we provide the full optimal weights for bispectrum analysis in the presence of sky-cuts and 
the inhomogeneous noise. The full inverse-covariance of the data is introduced which makes the estimator an optimal one. The final section 
is devoted to discussion and future plans for numerical implementation. 



2 MODELS OF PRIMORDIAL NON-GAUSSIANITY 

Deviations from pure Gaussian statistics can provide direct clues regarding inflationary dynamics. The single -field slow-roll model of inflation 
provi des a very small level of departure from Gaussianity, far below present experimental detection limits (Maldacena 2003; Acquaviva et al.l 
2003). Many other variants however will produce a much higher-level of non-Gaussianity which will be within the reach of all-sky experi- 
ments such as Planck. 

In general models can be distinguished by the way they predict coupling between different Newtonian potential modes: 
{^(kO^ka^ka)} = (27r) 3 5 3D (k! + k 2 + k 3 )-F(fci, fe, k 3 ). (2) 

The function F encodes the information about mode-mode coupling and 5 3D (ki + k 2 + k 3 ) ensures triangular equality in fc-space. Different 
models of inflation are typically divided into two different groups. The first class of model is known as the "local model", where the 
contribution from F{k\, fe, k 3 ) is largest when the wavevectors are in the so-called "squeezed" configuration, where e.g. k\ <C k 2 , k 3 . The 
local form of non-Gaussianity is predominant in models where there is non-linear coupling between the field driv ing inflation (the inflaton) 
and t he curvature perturbations (Salopek & Bo ng|l990lll99ll:lGangui et alj|l994|), such as the curvaton model ( iLvth. Ungarelli & Wandj 
120031) and the Ekpyrotic model ( Koyama e t alJl2007l :lBuchbinder . Khourv & Ovrut 2008). The primordial bispectrum in the local model can 
be written as: 



F loc (k 1 ,k 2 ,k 3 ) = f£i 



+ eye. perm. 



(3) 



where the power spectrum of inflationary curvature perturbations is given by Pq> — A$/fc 3 , in general for deviation from Harrison-Zeldovich 
power-spectra one has P$ = A* /fc 4- ™ 3 . 

The other main class consists of models where the contribution from F(ki,k 2 ,k 3 ) is maximum for configurations where ki ~ 
k 2 ~ ^3- The equilateral form appears from non-canonical kinetic terms such as Dirac-Born-Infield (DBI) action (Alishahiha et. al. 
2004), the ghost condensation (Akrani-Hamed e t al. 2004) or various single-field models where the scalar field acquires a low sound speed 
dChen. Easther & Liml2007l : ICheung et al.l20o3). The eq uilateral model is not a separable function of the ki, which complicates the analysis 



considerably, but it was shown by Creminel li et alj J2006t) that the following approximate form can model the equilateral case very accurately: 



F e "(k u k 2 ,k 3 ) = f^ L 



Al Al Al 

<pp - 2 T2pp + 6 T"T2T3 + eyeperm. 

1 2 1 2 3 2 3 



(4) 



Secondary non-Gaussianity resulting from various sources e.g. coupling of lensing and the Sunyaev-Ze l'dovich effect, or lensing and the In- 
tegrated Sachs-Wolfe effect can also contribute to the observed bispectrum ( Spergel & Goldberg 1999afi) . We will present a general analysis 
of sec ondary bispectra as well as the one induced by quasi-linear evolution of gravitational perturbations jMunshi. Souradeep & Starobinskvl 
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3 ANGULAR CMB BISPECTRUM WITH PRIMORDIAL NON-GAUSSIANITY 

The angular bispectrum can be defined as the three-point correlation in the harmonic domain. With temperature fluctuations as a function of 
solid angle fi, AT(fi), 

ai m = J dCl^^YiUfl) (5) 
and the three-point function may be written 

/ \-B ( I I' I" \ 

\0'lmO'l'm' a l"m" ) = &WI" „ I // ■ W 

\ m m m J 

This form preserves the the rotational invariance of the three-point correlation function in the harmonic domain. The quantity in parentheses 
is the Wigner 3j-symbol, which is nonzero only for triplets (Zi, h, I3) which satisfy the triangle rule, inclu ding that the sum h + h + h is 
even, ensuring the parity invariance of the bispectrum. The reduced bispectrum bun" was introduced bv lKomatsu & Spergell J200lh which 
will be helpful (see lBabich. Creminelli & Zaldarriaga] d2004l) for elaborate discussion): 



/ (2i + l)(2i' + l)(2i' + l) / I I' I" 

£>ll'l" = V : n n n I "If I" = Hl'l" "U'l" ■ U) 



4tt ^000 

The reduced bispectrum bun" can be expressed in terms of the kernel F(ki, k 2 , k^) for various models that we will be considering: 

b hhh = (^) 3 / drr2 J fcid*iiii(*ir)A^(fcir) J kldk 2 j l2 {k 2 r)Aj 2 {k 2 r) J fcidfcsJi^fcsrJA^CfcsrJFffci.fca.fcs); (8) 

where Af(k) denotes the transfer function which relates the inflationary potential $ to the spherical harmonics ai m of the temperature 
perturbation in the sky (e.g. IWang & Kamionkowskil fcOOOh '): 

i f d s k 



a lm = 4n(-i) 1 J — <£>(k)A/ (fe)F^(k). (9) 
Using these definitions one can express the reduced bispectra for the local and equilateral case as follows: 

b hhh = 2 fNL J r 2 dr [a h (r)f3i 2 (r)f3 h (r) + cyc.perm.] (10) 

b hhh = 6 fNL J r ' '' dr i~ a h ( r )Ph ( r ) A 3 (r) - 2S h (r)5 h (r)S h (r) + /3 h (r)-y h (r)S h (r) + cyc.perm.] (11) 

We will use these forms to construct associated fields A, B etc. from temperature fie lds with appropriate we ighting to optimise our estimator. 
We list the explicit expressions for the functions «;(r) , (3i(r) etc. for completeness dCreminelli et al. 2006): 



oti{r) = 


2 

7T 


/ fc 2 dfcAi(fc)ii(fcr) 
/o 


AM = 


2 

7T 


/ k 2 dkP<s,{k)Ai{k)ji{kr) 
Jo 


7;M = 


2 

7T 


f oo 

/ k 2 dkP# /3 (k)Ai(k)ji(kr) 
Jo 


Si(r) = 


2 

7T _ 


p oc 

/ fc 2 dfcP^ /3 (fc)A i (fc)j i (fcr) 
Jo 



(12) 

Numerical evaluations of these functions can be performed by using the publicly available software such as CAMB or CMBFAST. 



4 ALL SKY ANALYSIS WITH HOMOGENEOUS NOISE 

In this section, we compute the main statistic which will be used to estimate primordial non-gaussianity, and which can also be used to 
assess whether a non-gaussian signal is indeed primordial. We call the statistic the bispectrum-related power spectrum, as it derives from 
the cross-power spectrum of certain maps constructed from the CMB map data, and which, as we will see, is related to the primordial 
nongaussianity. 

The analysis in this section is optimal for detecting primordial non-gaussianity in the case of all-sky coverage and homogeneous noise. 
These assumptions will not hold in practice, but we present this simpler case for clarity first, and to show the connection with previous work. 
We relax the assumptions later, and give optimised estimators for realistic cases in later sections. 



4.1 Local Model 

Following Komatsu, Spergel & Wandeltl d2005h . we first construct the 3D fields A(r, O) and B(r, O) from the expansion coefficients of the 
observed CMB map, a; m . The harmonics here Ai m (r) and Bi m (r) are simply weighted spherical harmonics of the temperature field o; m 
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with weights constructed from the CMB power spectrum C'i and the functions a;(r) and Pi(r) respectively: 

A(r,h) = ^Y lm ((l)A lm (r)- Ai m {r) = ^-ha lm (13) 



/ rn 



Ci 



B(r,Cl) = y"Y lm (n)B lm (r); B lm (r) = ^P-b,a lm . (14) 

Im 

The function 6; represents beam sm oothing, and from here onward we will absorb it into the harmonic transforms. Using these definitions 
iKomatsu. Spergel & Wandeltl J2005^ define the one-point mixed-skewness for the fields A(r, Cl) and B(r, Cl): 

S l s oc = Sg B2 = J r 2 dr J dClA{r,Cl)B 2 {r,Cl). (15) 

5*3 can be used to estimate /jy£, but such radical data compression to a single number loses the ability to estimate contamination of the 
estimator by other sources of nongaussianity. As a consequence, we construct a less radical compression, to a function of I which can be used 
to estimate /]?£, but which can also be analysed for contamination by, for example, foregrounds. We do this by constructing the integrated 
cross-power spectrum of the maps A(r, Cl) and B 2 (r, Cl). Expanding B 2 in spherical harmonics gives 

B|m(*0 = / dClB 2 (r,Cl)Y lm (Cl) 



\- \- 01' (r) A" (r) 1 (21 + !)(«' + l)(2l" + 1) / I I' I" \ I I' I" \ 
1^2^ -QT^ry 4^ I m m' m" K™' a '"™" < 16 > 



and we define the cross-power spectrum Cf' B (r) at a radial distance r as 
C?' B '\r) = ^L J J2^{Aim(r)B^(r)] , (17) 

Integrating this over r gives: 

Cf B * = I r 2 drC?< B \r). (18) 



This integrated cross power spectrum of B 2 (r, Cl) and A(r, Cl) carries information about the underlying bispectrum -B;;';", as follows: 
C l' = ^TlEEE r dr \-cT Cv C V , \^a Vm ,a v , m „ 



m I'm' l"m n 



{ 2 l + mi > + m i» + i) i i> r w , v f _ 



47T \ / \ m m rn 

Similarly we can construct the cross power spectrum of the product map AB(r, Cl) and B(r, Cl), which we denote as C B ' AB ; 

c « = ^TlEEE J rdr {— Cl , c v , K a ""' a '"-" 

771 I'm' I" m" 



( 2 ; + i)( 2 , + i)(2^ + i) / z r r w , r r _ (20) 



47T VOOOyVmm m' 

Using these expressions, and the following relation, we can write this more compactly in terms of the estimated CMB bispectrum 

Bwi» = E ( m m' L"J ai ™ a ''™' a '"™" (21) 

mm' m" 

from which we compute our new statistic, the bispectrum-related power spectrum, C\ oc as 

n iloc , , f rjloc D 



or=vr +^ B ) = j0 T) EE{^^\ (22) 

where B\ v c v i is the bispectrum for the local /ml model, normalised to /j££ = 1. We can now use standard statis t ical t echniques to estimate 
/#£. Note that if we sum over all I values then we recover the estimator S pr i m of lKomatsu. Sperge l & Wandeh] J2005h . which is the cross- 
skewness of ABB: 

St = Si B2 =J2( 2l + ^ C t" 2 + 2 ^ S ' fl ) = ?hI E E E { B clcvCv' } ■ (23) 
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Figure 1. The bispectrum-related power spectrum is plotted as function of angular scale I. The red curve corrspond to the point source cross-contamination 
b pa = lO- 27 tiK and the bllue curve correspond to local model with /jvl = 1. See text for details. 

We show in Fig. 1 the form of the bispectrum-related power spectrum for the local model. If the signal observed is inconsistent with this form, 
it would indicate a departure from this form of primordial non-Gaussianity, and/or a significant contamination by foreground sources. In the 
right hand panel we show the expected form of contamination of the C\ oc statistic by a simple foreground source: unsubt racted point sources, 
rando mly distributed. The contamination scales with the number density. This contamination is expected to be very low (Komatsu & Spergel 
l200lh . but we show it as an illustration of how foreground effects may be detected. Since its structure is very different from the local 
primordial signal, it should be relatively easy to decouple. Note that the point source contamination here is the contribution to the local model 
bispectrum-related power spectrum - i.e. it is equation l !22t with one local B and one point source B (constant &;;'!")> not two point source 
B terms. 



4.2 Equilateral Model 

The form for the reduced bispectrum for the equilateral model as mentioned above is an approximation to the real bispectrum generated in 
theories with higher-order derivatives in the Lagrangian. In addition to the quantities A and B we have constructed analogous fields C and 
D with corresponding weights ji(r) and <5;(r): 




■ in 



(24) 



lm 




Im 



From these fields and using A and B previously defined a one-point statistic can be constructed dCreminelli et al. 2006): 




The associated power spectrum will have a composite structure with many contributing terms. 




Following the same procedure outlined above one can now write: 




(28) 



and finally we can recover the one-point statistic or the cross-skewness o f lKomatsu, Spergel & Wandeitl [2005). 




(29) 



Individual contributions to the final skewness following relations, which can be useful diagonistics for numerical checks: 





(30) 
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Given the signal-to-noise ratio of current estimates of S3 or equivalently /jvl from WMAP surveys, it may not be possible to use many 
narrow bins to evaluate the C;s associated with the primordial bispectra. However, with the increase in experimental sensitivity of future 
CMB experiments, it will be possible to divide the I range in narrower bins. The important point here is that the data must show consistency 
with the theoretical (j l ° c / eq j make a convincing case that the non-gaussianity is primordial. However, these expressions are only optimal 
for all-sky coverage and homogeneous noise; we relax these assumptions in the next sections. 



5 PARTIAL SKY COVERAGE AND NON-UNIFORM NOISE: AN APPROXIMATE TREATMENT 

It was pointed out in lBabichl d2005l) ; Icreminelli et al.l d2006l) : lYadav et al] d2008l) that in the presence of partial sky coverage, e.g. due to 
the presence of a mask or because of galactic foregrounds and bright point sources, as well as, in the case of non-uniform noise, spherical 
symmetry is destroyed. The estimator introduced above will then have to be modified by adding terms linear in the observed map. The linear 
terms for the local model can be written in the following form: 

S\Z ear = —^— I r 2 dr I dCl{2B(r,h){A(r,Cl)B(r,n)) sim +A(r,tl){B 2 (r,Cl)) sim } (31) 

fsky J J 

The linear terms therefore are constructed from correlating the Monte-Carlo (MC) averaged (A(n, r)B(n, r)) S i m product maps with the 
input B map. The mask and the noise that are used in constructing the Monte-Carlo averaged product map are exactly same as the observed 
maps and the ones derived from them such as A or B. 

S l : q near = -- f r 2 dr [ dQ{2B(r,Q){A(r,Q)B(n,r)) 3lrn + A(r,Cl){B 2 (r,n)) stm + 2D(r,Q){D(r,n) 2 ) 

fsky J J 

-2B(r, Cl) (C(r, Cl)D(Cl, r)) sim ~ 2C(r, fi) (B{r, Q)D(r, Cl)) sim - 2D(r, fi) (B{r, £l)C(r, &))}. (32) 

Mode-mode coupling is important at low angular modes, and we consider the full case later, but for higher frequency modes, we can 
approximate the linear correction to the local shape: 

C\° C = 1 \C A ' B2 - 2C ( - A ' B)B — C A '( B2 ' > X + 2 i(jAB,B _ p(AB),B _ £,B<A,B) _ gA(B,B) \ 
fsky v J fsky v J 

where f s k y is the sky fraction observed. 

The Cis such as C ; ( AS ^ S describe the cross-power spectra associated with Monte-Carlo averaged product maps {A(n,r)B(n,r)) 
constructed with the same mask and the noise model as the the observed map B. Likewise, the term qA(b,b) ( j enotes me avera g e cross- 
correlation computed from MC averaging, of the product map constructed from the observed map A(Q, r) multiplied with a MC realisation 
of map B(Q., r) against the same MC realisation B(Q., r). 



Ci q = -18 



±- {C^ 2 ~ 2C^ B - cf (B2) } + t~ {cf B B -C< 

Jsky ^ J J sky v 



(AB),B „B{A,B) n A(B,B) 



I -' / qD,D 2 _ 2(j(D,D)D _ £,D,<D 2 }~1 _ 2 (^B,CD _ 2(j(B,C)D _ ^B^CD)} 
" ky I 1 J fsky I ' J 



(34) 



2 ( „n.n 2 „^,(d.d)d ^D.tD 2 )^ 2 

fsky 

2 f fjD,CB _ 2(j(D,C)B gD.iCB) \ _ 2 f^,C,BD _ 2(j{C,B)D _ pC(BD) \ ' 
fsky v J fsky v J . 

ICreminelli etall ( 2006) showed via numerical analysis that the linear terms are less important in the equilateral case than in the local model. 
The use of such Monte-Carlo maps to model the effect of mask and noise greatly improves the speed compared to full bispectrum analysis. 

The use of linea r term s was found to greatly reduce the scatter of the estimator, thereby improving its optimality. The estimator was used 
in lYadav & Wandeltl |2008) also to compute the }nl from combined T and E maps. The analysis presented above for both the one-point 
statistics and the power- spectral analysis is approximate, because it uses a crude f s k y approximation to deconvolve the estimated power 
spectrum to compare with analytical prediction. A more accurate analysis should take into account the mode-mode coupling which can 
dominate at l ow I. The general expres sion which includes the mode-mode coupling will be presented in the next section. However it was 
found out by lYadav & Wandeltl d2008h that removing low Is from the analysis can be efficient way to bypass the mode-mode coupling. A 
complete numerical treatment for the case of two-point statistics such as Cf' B will be presented elsewhere. 



6 GENERALISATION OF OPTIMAL ESTIMATORS FOR REALISTIC SURVEY STRATEGY: EXACT ANALYSIS 

The general expression for the bispectrum estimator was developed bv lBabichl I200I) for arbitrary sky coverage and inhomogeneous noise. 
The estimator includes a cubic term, which by matched-filtering maximises the response for a specific type of input map bispectrum. The 
linear terms vanish in the absence of anisotropy but should be included for realistic noise to reduce the scatter in the estimates (see Babich 
2005 for details). We define the optimal estimator as: 

L' MM' ll'limm'mi ^ ' 

X { (^L'M',l 1 rn 1 a hrn 1 )(C lm l2m2 ai 2 rn 2 )(C l , m , l3m;j ai 3rn3 ) 
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' \m,V W '(fi L> M' ,l 2 m 2 a hm 2 ) 2C LMlm (^l'm',!jm2 a h™2 ) } 



(35) 



where N LL i is a normalization to be discussed later. A factor of 1/(21 + 1) can be introduced with the sum ^2 M , if we choose not to introduce 
the N LL i normalization constant. This will make the estimator equivalent to the one introduced in the previous section. Clearly as the data 
is weighted by C" 1 = (S + N) , or the inverse covariance matrix, addition of higher modes will reduce the variance of the estimator. In 
contrast, the performance of sub-optimal estimators can degrade with resolution, due to the presence of inhomogeneous noise or a galactic 
mask. However, a wrong noise covariance matrix can not only make the estimator sub-optimal but it will make the estimator biased too. The 
noise model will depend on the specific survey scan strategy. Numerical implementation of such inve rse- variance weighting or m ultiplication 
of a map by C -1 can be carried out by conjugate gradient inversion techniques. Taking clues from lSmifh & Zaldarriagal J2006I) . we extend 
their estimators for the case of the bispectrum-related power spectrum. We will be closely following their notation whenever possible. First 
we define Ql [a] and its derivative Ql [a] . The required input harmonics a; m are denoted as a. 

(L V I" \ 
M ml m" ) aVm ' ai " m " 
■ ... ., ... ' 



L I' I" \ >r- R L I I' 

I'm' ,1" m" M I'm' ^ 



(37) 



These expressions differ from that of one-point estimators by the absence of an extra summation index. Q l [a] therefore represents a map as 
well as di m QL[a], however Ql[<i] is cubic in input maps a; m where as di m QL[a] is quadratic in input. 

The bispectrum-related power spectrum can then be written as (summation convention for the next two equations): 

E L [a] = [N- 1 ]^ [QviC^a] - [C 1 ^(c^Q^C" 1 a']) mc)} (38) 

Here ()a/c denotes the Monte-Carlo averages. The inverse covariance matrix in harmonic domain l2Tri2 = (a; imi a; 2 , n ,) _1 encodes 

the effects of noise and the mask. For all sky and signal-only limit, it reduces to the usual C ( ~ mi i 2rrL2 = -^8ui8 mm i The normalisation of 
the estimator which ensures unit response can be written as: 

N LV = ~ ({d^Q^C-'a]} C^ miihm2 [d^Q^C^a]}) - i {(d^MC^a])} {{d^Q^C^a])} . (39) 
We will be using the following identity in our derivation: 

{[C _1 a]; imi [C _1 a]; 2m2 ) = C^ lii2m2 . (40) 

The Fisher matrix, encapsulating the errors and covariances on the El, for a general survey associated with a specific form of bispectrum 
can finally be written as: 

Fll ' = 12 12 B Ull ,B m ^ ( ^ ^ ^, ) ( ^' m 2 m> 2 

MM' l^'.m^m'. V ' V 

xi{2C" 1 f r, M ,C7 x ,, ,C7 X „ , +4C7, 1 , T , M ,C7 X , C7 X „ , ) = — (2aS + 4a*?, ) . (41) 

/? l_ L,M,L'M' l]_m\,l^m^ i2 m 2i'2 m 2 LM,L'M' ij_mi,t^m^ i2 rn 2 m 2 ^ ^ LL' J v J 

Using the following expressions which are extension of lSmith & Zaldarriagal J2006h . we find that the Fisher matrix can be written as sum of 
two a terms ce pp and a**. The alpha terms correspond to coupling only of modes that appear in different 3j symbols. Self couplings are 
represented by the beta terms. The subscripts describes the coupling of various I and L indices. The subscript PP correspond to coupling 
of free indices,i.e one free index L\ with another free index L 2 and similar coupling for indices that are summed over such as h, l 2 etc. 
Similarly for subscript QQ the free indices are coupled with summed indices. Couplings are represented by the inverse covariance matrices 
in harmonic domain e.g. CTs LM denotes coupling of mode LM with Im. 



Mi , hi 2 (ji'.mjm'. 



a L?L 2 ~ 12 12 B tihi[ B L 2 i 2 i> 2 f M \ m \ m '/ J f ]J 2 m 2 2 ^ J ClIm^l^C^^C^ ^ (42) 
12 B ^hi' 1 B L2i2i' 2 ( M \ mi mi ) ( M 2 m 2 m' 2 ) C ^Mi.h*» C hL 1 ,L a M a Ci' l m> 1 ji im ! l (43) 



a£f i2 = ( frjxpfojap ; a Q L % 2 = { Lihl^Ltl^ ) . (44) 

These results will reduce to those o f ISmith & Zaldarriagal d2006h when further summations over Li and L 2 are introduced to collapse the 
two-point object to the corresponding one-point quantity. The beta terms that denote cross-coupling can be written as: 



PLiL 2 - ^ B Llhlli B L2hl , 2 ^ ^ ^ ^, J 2 m 2 m' 2 ) C ^M 1 ,L 2 M 2 Cl 1 l ni ,l' 1 r n ' 1 C l 2 ln 2 ,l' 2 r n ' 2 ^ 



Mi,M 2 ! i !'.m i m< 
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Pl?l2 - Z Z B L 1 i 1 vB L ^ v _ i ( M ^ ^ m J , j f M ^ m 2 2 j C Xrl V 1 ,j 2m2 C'j i ^ ij j, mi C' Iia 1 Jf2> j, m , (46) 

Mi,M 2 lil'.mim'. v 7 v 7 

^il = £ £ ( Mi mi m\ ) ( Ma m 2 m 2 ) ^iWiim^Wia^jm^mJ (47) 

/3^ 2 = (LiWi)(L 2 j2ji); = ( £ 1 jij , i)(£ab 6); = (Liji p^j a) ■ (48) 

No summation over repeated indices is assumed. Using these expressions one can finally show that 

[F~ 1 }ll> = {E L E L> ) - (E L )(E L ,) = (AA + BB + CC + 2AB + 2BC + 2AC) LV (49) 
where 

AA LL , = { + ± a P[, + l/3£[ + ±pPP + ±prr } ; BB LL , = ; CC LL , = 4/3^, (50) 

2AB W = -2{pll + 202l)i 2A ^ll' = -4(20gg + 2BC 1L< = 4/3f« . (51) 

The final expression can be written in terms of only a terms as the j3 terms cancel out: 



2 



T? j " P P | ^ QQ 

1 36 LL 36 M 



'}• (52) 

If we sum over LL' the Fisher matrix reduces to a scalar F = Y\ v Fll' w ith, = <*2l' = Q an ^ /^lT' = /^fx/ ~ ^li?' ~ P' wnere 
a, /3 and F are exactly the same as introduced in ISmith & Zald arriaga (2006). 



6.1 Joint Estimation of Multiple bispectrum-related Power-Spectra 

The estimation technique described above can be generalised to take cover the bispectrum-related power spectrum associated with different 
set of bispectra (X,Y): 

E L [a] = [F- 1 ]^, {QUC^a} - [C' 1 a] lm {d lm Q Y L , [C^a^Mc) } . (53) 

The associated Fisher matrix now will consist of sectors F£~* ,F Y L , and F LL , , The sector XX and YY will in general will be related to 
errors associated with estimation of bispectra of X and Y types, whereas the sector XY will correspond to their cross-correlation. 

r If f 2 r ppiXY 4 qq, X y\ r<AS 

LL ' = 136 ^ LL '\ 36 [ l v\ } ( 54 ) 

where we have: 

= Z B *hi[ B L>i 2 i'J m ^ m'j j ( M' m 2 m 2 ) Cf iAf,£'M' C U^,l a m a ^m' I ,JJmJ ^ 

MM' lil'.mim'. V 7 V 7 

and a similar expression holds for [oi^,] XY ■ 



6.2 Generalisation to non-optimal weights 

Although the estimator as constructed is fully optimal - its true usefulness is determined by the affordability of the construction of the C _1 
matrix as well as the availability of a fast method to multiply it with CMB maps in a Monte Carlo chain. A more general class of estimator 
which is sub-optimal can be constructed by replacing the inverse covariance weighting of the data [C~ a] by [Ra], where [R] is an arbitrary 
filter function. In this case the estimator with unit response can be written as: 

E?[a] = [^ll' [QL'\Ra] - [Ra]i m (d lm Q L >[Ra\) MC ) (56) 

L' 

where the normalisation and the variance of the estimator can be constructed in a similar manner: 

F LL ' = ({EL^E^-iE^iEv)- 1 ^^({dtmQLlRsVRidtmQL'lRs^MC 

= \({di imi Qi[Ra]}[FCF] limi ,i 2m2 {d limi Qv[Ra]}) - | {{d^Q^Ra])} C^M^ {{d h m 2 Qi>[Ra])} . (57) 

The optimal weighting can be replaced by an arbitrary weight or no weighting at all (R=I). However in this case the estimator though unbiased 
clearly becomes a sub-optimal one. 
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6.3 Recovery of all-sky homogeneous noise model 

In the the all-sky limit we recover the usual expression: 

I'll- = ^ <1 > ' 5 LL , +4V \ . (58) 



36 I ^ di 

\ iv 



36 ^ CiClCv CiClC l , 
In the case of joint analysis as before we can write down the off-diagonal blocks of the Fisher matrix as: 



p XY _ 1 J o B LlV B Lll> c | a \ " B LL'l B LL'l ( . 

F "'- 36 2 L CMC, S ^'+ 4 2^ Cl C l ,C x > - m 
I w i 



For X — Y we recover the diagonal blocks of the Fisher matrix for the independent estimations derived before. The errors for independent 
estimates are given by */ {F**)- 1 , where as for the joint estimation the errors are U (F\ 



LL I 



6.4 More general bispectra 

A more general bispectrum can be written as a sum of individual product terms: 

N fact 

6 



1 iact 

bhhh = 5 A h B h C h +symm.perm. (60) 



This can be seen as a generalization of the type of bispectra introduced for the equilateral case. It may also possibl e to approximate the 
bispec trum bi 1 i 2 i 3 to a smaller number of optimum factorizable terms N opt with suitable weight factors Wi as given in lSmith & Zalda rriaga 
(2006). In this case the bispectrum is expressed as: 



N„ 



b hhh = g X! WiA h B h C h + symm.perm. (61) 



Clearly significant computational gain can only be achieved if N op t "C Nf ac f In the following discussion we generalise the description in 
previous sections to such composite bispectra. Following the same analytical reasoning we can show that the Fisher matrix elements for such 
a composite bispectrum can be written as: 

(2L + l)(2/ + l)(2/' + l) ( L I I 1 N 2 



1f xy } JL y (2L + l)(2i + l)(2r + l) / L I V \ 1 r t j j ,xr ^ , f 

\F LL ,] l3 - -\25 LL , — I Q Q Q I ^ d ^[A L B l C l , + ...\ [[A^C. + .-.l 

iv v 7 



(62) 

The symbols A 1 , B % etc denotes the i-th term in the factorised representation of a specific type of the bispectrum of type X or Y. The total 
contribution of all terms will constitute the final Fisher matrix: 

F% = $>fj] y ; F XY = Y, (63) 

ij LL' 

Using the following identity we can project this expression onto real-space: 

dzP h (z)P l2 (z)P h (z) = 2f jj h Q M (64) 



(65) 



iFLL'h =j *> [(& A3 { Z )e BJ ( Z )e cJ (z) + . . .) s LLI + (tf Ai ( Z )d: Bi (*)e cl w + . . .) 

The first term describes the diagonal entries of the Fisher matrix and the second term relates to the off-diagonal terms. 

i 

In case of cross-correlational studies the A 1 and A 3 will come from two different factorization of distinct bispectra denoted before by X and 
Y. The case of weighted sum can also be derived in an exactly similar manner. 



7 MODELS FOR NON-GAUSSIANITY 

We will use two commonly-used specific models for the primordial non-Gaussianity as well as one foreground source of contamination i.e. 
extra-galactic point sources to demonstrate the power of our statistics in this section. 
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7.1 Local or squeezed Model 

Using the specific form for ftj^jg the Fisher matrix elements for the local model can be expressed in terms of the functions a and (3 and the 
power spectrum Ci . 

2 



a 



o o o ; c Ll c La Ci 

y { I r',lr(2n L ( r) i Lj (r) l,(r) + ai(r)/3 Ll (r)f3 L2 (r))\ (67) 



Zi : - & lia ^(2Lr + 1) J2(2h + 1)(2I 2 + 1) ( Ll h ' 2 " " 



1L2 ^ 11 a«i n „ o J C Ll C h d 



|y r 2 dr (a Ll {r)p h (r)/3, a (r) + a h (r)A 2 (r)^ 1 (r) + a, 2 (r) (3 h {r)(3 Ll (r)) | . (68) 



7.2 Equilateral Model 

Using the equilateral model the contribution to the Fisher matrix can be expressed as: 



<l 2 = %(2Lr + 1)(2L 2 + 1) £(2J + 1) ( ^ h Q 

1 ^ 



t2 



2 



1 



/ Cl 1 Cl 2 Ci 



r 2 dr (-a Ll (r)/3 L2 (r)f3i(r) + 8 Ll (r)8 L2 (r)8i(r) + f3 Ll (r)j L2 (r)8i(r) + cyc.perm.) > . (69) 



^ 2 ^ 1 + i)^E( 2 ^ 1 )( 2 ^ 1 )( L o 1 h o) 2 cI^c l 



r 2 dr (-a Ll (r)/3 h (r)/3 h (r) + S Ll (r)S h (r)8 h (r) + (3 Ll (r) 7il (r) 7i2 (r) + cyc.pcrm.) | (70) 

The power spectra Ci appearing in the denominator take contributions both from the pure signal or CMB and the detector noise. It is possible 
to bin the estimates in large enough bins to report uncorrected estimates which may be possible for an experiment such as Planck with very 
high sky-coverage. A detailed analysis of the singularity structure of the error-covariance matrix will be presented elsewhere. 



7.3 Point Sources 

The bispectrum from residual point-sources which are assumed random can be modelled as bi x i 2 i 3 — b ps . The exact value depends on the 
flux limit and the mask used in the survey. The accuracy of such an approximation can indeed be extended by adding contributions from 
correlation terms. 



pp 
a Ll L 2 



a 7 T 



Similar computations can be done for cross-correlation among various contributions, e.g. contamination due to point sources or estimation 
of a specific type of non-Gaussianity. A joint estimation is useful for finding out also the level of cross-contamination from one theoretical 
model while another is being estimated. 



8 CONCLUSIONS 

We have addressed the problem of finding an estimator of primordial non-Gaussianity from microwave background data. The new feature 
of this analysis is that the technique presented here allows one to make an assessment of whether any non-Gaussian signal is primordial or 
not. The issue here is that if one finds an estimate of a level of primordial non-Gaussianity which is inconsistent with zero, then it is very 
difficult at present to make a convincing case that it is indeed primordial and not simply contamination by any number of other effects which 
might lead to a non-Gaussian CMB map. The method does this by performing a less aggressive data compression than previous analyses. 
Rather than compressing the data to a single number (typically an estimate of /jvl), it reduces the data to a function, the bispectrum-related 
power spectrum. This is an average cross-power spectrum of certain maps constructed from the CMB data. By doing this construction, 
one retains the ability to assess the contributions from different sources, such as residual point sources, incomplete foreground subtraction 
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and so on (see e.g. lSerra & Cooravl J2OO8)). As an example, we have computed the expected bispectrum-related power spectrum optimised 
for local non-gaussianity (we also consider the equilateral type), and calculated the contribution expected from an unclustered population 
of point sources. Indeed, one can use standard statistical methods to estimate the amplitude of components of non-Gaussianity, and since 
these contributors have quite different harmonic dependences, the estimators will be largely decoupled. T he power of the technique will 
depend on the level of the primordial signal, but if it is at the level claimed by Yadav & Wandelt ( 2008), then it will be possible with 
Planck data to construct a large number of band-power estimates of the bispectrum-related power spectrum to see if it is primordial. We 
also include polarisation in the ana lysis (see appendix). The work draws extensively on previous studies, in particular generalising the 
iKomatsu. Spergel & WandeldJiibll) analysis for the all-sky, homogeneous noise case. For the more realistic case of partial sky coverage and 
inhomogeneous noise, we present optimised estima t ors of the bispectrum-related power sp ectrum, including linear terms and extending the 
work of lBabichl d2005l) : ISmifh & Zaldarriagi] j2006h ; ISmithT Scnatorc & Zaldarriagi] f2009). For studies such as this, it is normal to assume 
the b ackground cosmology is known from the power spectrum, but uncertainties will propagate into the Jnl estimates dLigouri & Riottd 
2008) and will be investigated elsewhere. Note that the techniques here could be generalised to higher-order statistics such as the trispectrum, 
should the bispectrum vanish for some symmetry reason. 
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APPENDIX A: JOINT ANALYSIS OF TEMPERATURE AND POLARIZATION 

Most current constraints on non-Gaussianity still come from temperature maps, but with WMAP and Planck, the situation is changing. 
Similar calculations can be performed for joint temperature and E-type polarisation analysis, and the estimators discussed above can be 
generalised to include E -type polarisation to tighten the constraints. The functions on and f3 i that we discussed in the main tex t needs to 
be generalised for both T, t e mperature and i?-ty pe polarization. We follow the discussion in lYadav, Komatsu.& Wandeltl d2007h . see also 
Babic h & Zaldarriaga! d2004l) ; lLiguori et all j2007l) for related discussions. 

2 f°° 2 f°° 

af(r) = - k 2 dkP^,(k)Af(k)ji(kr) (r) = -/ k 2 dkA? (k)j,(kr). (Al) 

* Jo n Jo 

The index X can be T or E. Similar calculations can in principle be performed for the equilateral case. We will focus however here only 
on the local or squeezed model. The power spectra now can be a temperature (TT) -only power spectrum or an Electric-Electric (EE) 
polarization spectra which we define below, with specific choice of A;. 

2 f°° 

C? Y (r) = - / k 2 dkP<s,(k)A?(k)AY(k). (AT) 
71 Jo 

We arrange the all-sky covariance matrix (in the harmonic domain) in a matrix, which takes the following form in terms of the Ci defined 
above: 



[C\i 



(jTE (J EE 



(A3) 



We can now construct the A(r, ft) and B(r, ft) fields now can be constructed using the full covariance matrix instead of only temperature 
data. 



A(r, ft) = E 7 <*'ma?(r)Yjm(A) (A4) 

lm ip 

B(r,Cl) =J2JL [ C "T «<m#(r)y,m(ft). (A5) 



lm ip 



The construction now follows exactly the same steps as depicted for the temperature case. We compute the cross-correlation of B 2 (r, ft) 
with the A(r, ft). The B 2 (r, ft) can be decomposed in terms of the T and E harmonics both a\ m (j here takes values T or E). 

CW= / dnB 2 (r,£l)Y lm (n) 



ST \T a r\ ( \\n-^ kT / (2Z + 1)(2Z^ + 1)(2Z" + 1) f I I' I" \ ( I I' I" \ 3 k 

I'm' l"m" V ' V ' 

The cross-correlation of the associated A and B field will contain information both from temperature and polarization maps. 

C?' b2 = ( r 2 drC^ B \r) = I r 2 dr Bf m (r) A lm (r) = 



J r 2 drBf m (r)A lm (r) = ^ E E E / [ C ~T ^ [^V ^ [ C ~V' 

m I'm' I" m" 



(2/ + l)(2Z' + l)(2/"+l) / I I' I" \ ( I I' I 



,_ 1 (, „ o / \ m m' m" ) a> ^ a i'm' a i"m" (A7) 

i:;/: r , = [ l m l m > !n" ) a ' "i',„'"'i" " (A8i 

mm' m" 
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The mixed-bispectrum Bf.JL, , contains information about three-point correlation in harmonic space with the possibility that the harmonics 
a\ m can either be of T or temperature type or E or electric-polarisation type. The theoretical model for such bispectrum depends on functions 
af (r) and f3 q (r) where again depending on the superscript the functions can be of temperature or the electric type. 



Bfiu" = ^ W + lW + V^ l+li ( J J ^ J r * dr {ft {r )ft, (r)aF« (r) + # (r)a«, (r)tf„ (r) + of (r)#, (r)tf„ (r) } .(A9) 
Finally in an analogous way we can express the power- spectrum related to the mixed-bispectrum as follows: 

(2/ + l)(Cr 2 +2C? B > B ) = F L Y,Y,{ B ^ [C-T [C- 1 ]'; [C^B^,,}. (MO) 



V V' 



We have assumed summation of repeated indices which denote the polarisation types, i.e. k and p, q, r in the above expression. The 
one-point mixed skewness which is related to the above power spectrum is analogously written as follows: 

5--=^(2 i + l)(C 2 + 2C B )=f^^^{^ [C-X [C-'Y; [C-'t Btf„} . (All) 



i v 



This generalises the temperature-only power spectrum estimator introduced in the text of this paper and can be extended to take into account 
other models, sky-coverage and secondary anisotropies in an analogous manner. 



